Jack Colvin, Along with my group, Larrence Xing, Kevin Wu, Annabelle Yarborough and Olivia Lucal
Sersic Profiles, Introduced by Argentinian astronomer José Luis Sérsic, are a generalization of de Vacouleurs' law, which describes the light prodiles of elliptical galaxies, to all galaxy types. It specically fits galaxies to a symmetrical light distribution described by the equation $I * e^{-b((\frac{r}{r_{e}})^{\frac{1}{n}}-1)}$, where n is a varible index that describes how quickly the light in a given galaxy distribution falls off from its max value, with a higher index representing a sharper fall-off. This is useful in studying galaxies specifically in how spiral and elliptical galaxies are sharply dilineated with respect to this index: eliptical galaxies generally have sersic index of around 4, while spiral galaxies have an index of around 1 (which is an exponential profile). With this in mind, I took 6 observations, 3 of spiral galaxies (m51, m81, m101) and three eliptical gal;axies (m89, m49, m60) to see if I could
I started by taking 120s exposures of each galaxy with seo, which represented a good balance betweenb having a high enough SNR, while avoiding the pixels within the galaxies from spilling into neigboring pixels (which would cause big problems whgen trying to conduct sersic fits). I took exposures in all three bands in order to be able to create pretty rgb images for all three galaxies! Finaly, I constructed a sersic function and a least sqaures residual function in order to use least sqaures fitting to fit a sersic model to each galaxy.
While this relatively simple model worked for elliptical galaxies relatively well, I had to make some modifications when working with spiral galaxies. First and foremost, the bulge is both much brighter than the rest of the disk, and has a different, more concentrated sersic profile, which was resulting in a one-component model only fitting to the bulge, and not to the rest of the galaxy, and thus overestimating the overall sersic index. To correct for this, I switched to a two-component model when dealing with disks. I also suspect that the spiral arms are also caussing fitting issues, as they are assymetric, and are likely thus ignored by the model. I hav3 enot yet corrected for this effect, but I will certainly do so as I continue to work on this project in the future. Finally, I tried masking out the bulge entirely, to see if doiung so would give me a better or worse result than the two component model, with varying results, and I also tried fitting a two component model with a sersic index for the bulge set between 3-6 as another method to remove the variability of the bulge as my goal is to fit a model to the disk halo.
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.visualization import ZScaleInterval
from matplotlib.colors import LogNorm
from astropy.nddata import Cutout2D
from astropy.modeling import models, fitting
from astropy.visualization.lupton_rgb import make_lupton_rgb
from astropy.visualization import LogStretch
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
def plot_prettier(dpi=150, fontsize=11, usetex=False):
'''
Make plots look nicer compared to Matplotlib defaults
Parameters:
dpi - int, "dots per inch" - controls resolution of PNG images that are produced
by Matplotlib
fontsize - int, font size to use overall
usetex - bool, whether to use LaTeX to render fonds of axes labels
use False if you don't have LaTeX installed on your system
'''
plt.rcParams['figure.dpi']= dpi
plt.rc("savefig", dpi=dpi)
plt.rc('font', size=fontsize)
plt.rc('xtick', direction='in')
plt.rc('ytick', direction='in')
plt.rc('xtick.major', pad=5)
plt.rc('xtick.minor', pad=5)
plt.rc('ytick.major', pad=5)
plt.rc('ytick.minor', pad=5)
plt.rc('lines', dotted_pattern = [2., 2.])
if usetex:
plt.rc('text', usetex=usetex)
else:
plt.rcParams['mathtext.fontset'] = 'cm'
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman'] + plt.rcParams['font.serif']
plot_prettier(dpi=150, fontsize=11)
iband = fits.open("/Users/jackcolvin//m51_i-band_120.0s_bin1_250603_054315_jackcolvin_seo_0_FCAL.fits")
gband = fits.open("/Users/jackcolvin//m51_g-band_120.0s_bin1_250603_054042_jackcolvin_seo_0_FCAL.fits")
rband = fits.open("/Users/jackcolvin//m51_r-band_120.0s_bin1_250603_053809_jackcolvin_seo_0_FCAL.fits")
r = iband[0].data
g = rband[0].data
b = gband[0].data
#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled = zscale(r)
g_scaled = zscale(g)
b_scaled = zscale(b)
image = make_lupton_rgb(r_scaled, g_scaled, b_scaled, Q = 1, stretch = .9, minimum = .3)
plt.imshow(image, origin = 'lower')
plt.show()
#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction
#based on whether this galaxy is spiral or elliptical.
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r, (1100, 1200), (1000, 1000))
data0 = cutout.data
data = (data0-np.median(r))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )
y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
#define an example sersic profile:
def sersic_2d(x, y, I_e, R_e, n, x0, y0):
r = np.sqrt((x - x0)**2 + (y - y0)**2)
b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
return I_e * np.exp(-b_n * ((r / R_e)**(1/n) - 1))
def residuals(params, x, y, data):
I_e, R_e, n, x0, y0 = params
model = sersic_2d(x, y, I_e, R_e, n, x0, y0)
return (model - data).ravel()
p0 = [np.max(data), 50, 1.0, 500, 500] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")
Fit: I_e = 0.04911, R_e = 194.49, n = 2.33, x0 = 493.56, y0 = 500.00
model_image = sersic_2d(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax )
plt.show()
#we can now evaluate the chi2 of our fit:
data1 = data.flatten()
mask = data1 > 0
data2 = data1[mask]
model_image1 = model_image.flatten()
model_image2 = model_image1[mask]
residual = (data2 - model_image2)
uncertainty = (np.sqrt(data2))#use poisson uncertainty on fluxes
chi2 = np.sum((residual/uncertainty)**2)
dof = len(data2) - len(result.x)
chi2_red = chi2/ dof
print(f"Chi-squared: {chi2:.2f}")
print(f"Reduced Chi-squared: {chi2_red:.2f}")
#very low reduced chi2, almost symptomatic of over fitting
#also, our sersic index is very high, almost 4
#which indicates that it is capturing only the bulge, rather than the full disk
#so now I wanted to try a two-component model
Chi-squared: 59764.06 Reduced Chi-squared: 0.10
def two_component_model(x, y, params):
# params = [I_e_bulge, R_e_bulge, n_bulge, I_e_disk, R_e_disk, x0, y0]
I_e_bulge, R_e_bulge, n_bulge, I_e_disk, R_e_disk, n_disk, x0, y0 = params
bulge = sersic_2d(x, y, I_e_bulge, R_e_bulge, n_bulge, x0, y0)
disk = sersic_2d(x, y, I_e_disk, R_e_disk, n_disk, x0, y0)
return bulge + disk
def residuals_two_components(params, x, y, data):
model = two_component_model(x, y, params)
return (model - data).ravel()
p1 = [np.max(data), 50, 4, np.max(data)/1.5, 200, 1, 500, 500]
bounds_lower = [0, 1, .5, 0, 10, .5, 0, 0]
bounds_upper = [np.inf, 200, 8, np.inf, 500, 8, data.shape[1], data.shape[0]]
result = least_squares(residuals_two_components, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))
I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, x0_fit, y0_fit = result.x
print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")
Fit results: Bulge: I_e=0.72, R_e=18.74, n=1.05 Disk: I_e=0.07, R_e=185.16, n=0.77 Center: x0=494.71, y0=515.74
model_image = two_component_model(x, y, np.array(result.x))
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#thus masking pout the spiral arms may be a necessary step to getting a better fit.
#now try fitting to the disk profile by removing the center bulge entirely
#to see if the fit becomes becomes better after doing so
from astropy.convolution import Gaussian2DKernel, convolve
kernel = Gaussian2DKernel(x_stddev=20) # tune stddev to smooth scale
smooth = convolve(data, kernel)
residual = data - smooth
sigma = np.std(residual)
threshold = sigma
R_bulge = 50 # adjust based on your galaxy
yy, xx = np.indices(data.shape)
r = np.sqrt((xx - 500)**2 + (yy - 500)**2)
central_mask = r < R_bulge
# Mask = True where pixels are *not* spiral arms (we want to keep them)
#but we also need to avoid masking out the central bulge.
mask0 = np.abs(residual) < threshold
mask = mask0
data_corr = np.where(mask, data, np.median(data))
y = np.zeros(shape = data_corr.shape, dtype = int)
x = np.zeros(shape = data_corr.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
p1 = [np.max(data), 50, 4, np.max(data)/1.5, 200, 1, 500, 500]
bounds_lower = [0, 1, .5, 0, 10, .5, 0, 0]
bounds_upper = [np.inf, 200, 8, np.inf, 500, 8, data.shape[1], data.shape[0]]
result = least_squares(residuals_two_components, p1, args=(x, y, data_corr), bounds=(bounds_lower, bounds_upper))
I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, x0_fit, y0_fit = result.x
print(f"Smoothed Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")
Smoothed Fit results: Bulge: I_e=0.00, R_e=1.00, n=0.50 Disk: I_e=0.06, R_e=178.04, n=0.84 Center: x0=486.89, y0=529.90
model_image = two_component_model(x, y, np.array(result.x))
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Smooth_data')
plt.imshow(data_corr, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data_corr - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#While the bulge plot is no longer accurate
#as the bulge is filtered out
#the model sersic fit is much closer to what we expect
#and the residual plot looks much bteter for the outside.
iband.close()
rband.close()
gband.close()
iband0 = fits.open("/Users/jackcolvin//m101_i-band_120.0s_bin1_250603_053348_jackcolvin_seo_0_FCAL.fits")
gband0 = fits.open("/Users/jackcolvin/m101_g-band_120.0s_bin1_250603_053050_jackcolvin_seo_0_FCAL.fits")
rband0 = fits.open("/Users/jackcolvin//m101_r-band_120.0s_bin1_250603_052816_jackcolvin_seo_0_FCAL.fits")
r0 = iband0[0].data
g0 = rband0[0].data
b0 = gband0[0].data
#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled0 = zscale(r0)
g_scaled0 = zscale(g0)
b_scaled0 = zscale(b0)
image = make_lupton_rgb(r_scaled0, .9*g_scaled0, b_scaled0, Q = 1, stretch = .7, minimum = .25)
plt.imshow(image, origin = 'lower')
plt.show()
#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction
#based on whether this galaxy is spiral or elliptical.
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r0, (1100, 1250), (600, 600))
data0 = cutout.data
data = (data0-np.median(r0))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )
y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
#define an example sersic profile:
def sersic_2d(x, y, I_e, R_e, n, x0, y0):
r = np.sqrt((x - x0)**2 + (y - y0)**2)
b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
return I_e * np.exp(-b_n * ((r / R_e)**(1/n) - 1))
def residuals(params, x, y, data):
I_e, R_e, n, x0, y0 = params
model = sersic_2d(x, y, I_e, R_e, n, x0, y0)
return (model - data).ravel()
p0 = [np.max(data), 50, 1.0, 250, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")
Fit: I_e = 0.00152, R_e = 2192.04, n = 4.42, x0 = 237.36, y0 = 298.53
model_image = sersic_2d(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax )
plt.show()
p1 = [np.max(data), 10, 4, np.max(data)/1.5, 200, 1, 300, 280]
bounds_lower = [0, 1, 2, 0, 10, .5, 0, 0]
bounds_upper = [np.inf, 200, 8, np.inf, 500, 8, data.shape[1], data.shape[0]]
result = least_squares(residuals_two_components, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))
I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, x0_fit, y0_fit = result.x
print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")
Fit results: Bulge: I_e=0.06, R_e=48.55, n=2.00 Disk: I_e=0.04, R_e=270.66, n=0.92 Center: x0=237.33, y0=298.52
model_image = two_component_model(x, y, np.array(result.x))
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#thus masking pout the spiral arms may be a necessary step to getting a better fit.
def two_component_ellipse(x, y, params):
# params = [I_e_bulge, R_e_bulge, n_bulge, q_bulge, theta_bulge, I_e_disk, R_e_disk, q_disk, theta_disk, x0, y0]
I_e_bulge, R_e_bulge, n_bulge, q_bulge, theta_bulge, I_e_disk, R_e_disk, n_disk, q_disk, theta_disk, x0, y0 = params
bulge = sersic_2d_elliptical(x, y, I_e_bulge, R_e_bulge, n_bulge, x0, y0, q_bulge, theta_bulge)
disk = sersic_2d_elliptical(x, y, I_e_disk, R_e_disk, n_disk, x0, y0, q_disk, theta_disk)
return bulge + disk
def residuals_two_ellipse(params, x, y, data):
model = two_component_ellipse(x, y, params)
return (model - data).ravel()
p1 = [np.max(data), 15, 4, .5, 1, np.max(data)/2, 200, 1, .5, 1, 230, 270]
bounds_lower = [0, 0, 2, 0, 0, 0, 50, .5, 0, 0, 0, 0]
bounds_upper = [np.inf, 50, 6, 1, 2*np.pi, np.inf, 400, 8, 1, 2*np.pi, data.shape[1], data.shape[0]]
result = least_squares(residuals_two_ellipse, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))
I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, q_bulge_fit, theta_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, q_disc_fit, theta_disc_fit, x0_fit, y0_fit = result.x
print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}, q ={q_bulge_fit:.2f}, theta = {theta_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}, q={q_disc_fit:.2f}, theta={theta_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")
Fit results: Bulge: I_e=0.07, R_e=50.00, n=2.00, q =0.81, theta = 5.98 Disk: I_e=0.03, R_e=313.05, n=1.01, q=0.80, theta=2.54 Center: x0=237.30, y0=298.50
model_image = two_component_ellipse(x, y, np.array(result.x))
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#thus masking pout the spiral arms may be a necessary step to getting a better fit.
iband0.close()
gband0.close()
rband0.close()
iband01 = fits.open("/Users/jackcolvin//m81_i-band_120.0s_bin1_250603_055612_jackcolvin_seo_0_FCAL.fits")
gband01 = fits.open("/Users/jackcolvin//m81_g-band_120.0s_bin1_250603_055340_jackcolvin_seo_0_FCAL.fits")
rband01 = fits.open("/Users/jackcolvin//m81_r-band_120.0s_bin1_250603_055109_jackcolvin_seo_0_FCAL.fits")
r01 = iband01[0].data
g01 = rband01[0].data
b01 = gband01[0].data
#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled01 = zscale(r01)
g_scaled01 = zscale(g01)
b_scaled01 = zscale(b01)
image = make_lupton_rgb(r_scaled01, .85*g_scaled01, b_scaled01, Q = 1, stretch = .7, minimum = .25)
plt.imshow(image, origin = 'lower')
plt.show()
scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r01, (1000, 1200), (1000, 1000))
data0 = cutout.data
data = (data0-np.median(r01))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )
y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
#define an example sersic profile:
def sersic_2d(x, y, I_e, R_e, n, x0, y0):
r = np.sqrt((x - x0)**2 + (y - y0)**2)
b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
return I_e * np.exp(-b_n * ((r / R_e)**(1/n) - 1))
def residuals(params, x, y, data):
I_e, R_e, n, x0, y0 = params
model = sersic_2d(x, y, I_e, R_e, n, x0, y0)
return (model - data).ravel()
p0 = [np.max(data), 50, 1.0, 500, 500] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 600, 600]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")
Fit: I_e = 0.20730, R_e = 148.97, n = 3.14, x0 = 501.55, y0 = 491.44
model_image = sersic_2d(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax )
plt.show()
#a circularly symetric model does a bad job for this galaxy
#so now we must define an ellipse
def sersic_2d_elliptical(x, y, I_e, R_e, n, x0, y0, q, theta):
x_shift = x - x0
y_shift = y - y0
x_rot = x_shift * np.cos(theta) + y_shift * np.sin(theta)
y_rot = -x_shift * np.sin(theta) + y_shift * np.cos(theta)
R = np.sqrt(x_rot**2 + (y_rot / q)**2)
b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
return I_e * np.exp(-b_n * ((R / R_e)**(1/n) - 1))
def residuals_elliptical(params, x, y, data):
I_e, R_e, n, x0, y0, q, theta = params
model = sersic_2d_elliptical(x, y, I_e, R_e, n, x0, y0, q, theta)
return (model - data).ravel()
p0 = [np.max(data), 50, 1.0, 500, 500, .5, np.pi/3] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals_elliptical, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0, 0, 0], [np.inf, np.inf, 10, 600, 600, 1, 2*np.pi]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit, q_fit, theta_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}, q = {q_fit:.2f}, theta = {theta_fit:.2f}")
Fit: I_e = 0.20369, R_e = 181.48, n = 3.14, x0 = 501.58, y0 = 491.36, q = 0.71, theta = 0.99
model_image = sersic_2d_elliptical(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax )
plt.show()
p1 = [np.max(data), 50, 4, .5, 1, np.max(data)/1.5, 200, 1, .5, 1, 500, 500]
bounds_lower = [0, 1, .5, 0, 0, 0, 10, .5, 0, 0, 0, 0]
bounds_upper = [np.inf, 200, 8, 1, 2*np.pi, np.inf, 500, 8, 1, 2*np.pi, data.shape[1], data.shape[0]]
result = least_squares(residuals_two_ellipse, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))
I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, q_bulge_fit, theta_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, q_disc_fit, theta_disc_fit, x0_fit, y0_fit = result.x
print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}, q ={q_bulge_fit:.2f}, theta = {theta_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}, q={q_disc_fit:.2f}, theta={theta_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")
Fit results: Bulge: I_e=3.00, R_e=14.68, n=0.98, q =0.83, theta = 0.71 Disk: I_e=0.22, R_e=206.61, n=1.75, q=0.60, theta=1.06 Center: x0=501.58, y0=491.31
model_image = two_component_ellipse(x, y, np.array(result.x))
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
p1 = [np.max(data), 50, 5, .5, 1, np.max(data)/2, 200, 1, .5, 1, 500, 500]
bounds_lower = [0, 10, 2.5, 0, 0, 0, 0, .5, 0, 0, 0, 0]
bounds_upper = [np.inf, 50, 6, 1, 2*np.pi, np.inf, 500, 8, 1, 2*np.pi, data.shape[1], data.shape[0]]
result = least_squares(residuals_two_ellipse, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))
I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, q_bulge_fit, theta_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, q_disc_fit, theta_disc_fit, x0_fit, y0_fit = result.x
print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}, q ={q_bulge_fit:.2f}, theta = {theta_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}, q={q_disc_fit:.2f}, theta={theta_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")
Fit results: Bulge: I_e=1.33, R_e=10.00, n=2.50, q =0.02, theta = 3.22 Disk: I_e=0.20, R_e=182.42, n=3.14, q=0.71, theta=1.00 Center: x0=501.57, y0=491.45
model_image = two_component_ellipse(x, y, np.array(result.x))
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
iband01.close()
gband01.close()
rband01.close()
iband1 = fits.open("/Users/jackcolvin//m89_i-band_120.0s_bin1_250603_045958_jackcolvin_seo_0_FCAL.fits")
gband1 = fits.open("/Users/jackcolvin//m89_g-band_120.0s_bin1_250603_045721_jackcolvin_seo_0_FCAL.fits")
rband1 = fits.open("/Users/jackcolvin//m89_r-band_120.0s_bin1_250603_045449_jackcolvin_seo_0_FCAL.fits")
from astropy.visualization.lupton_rgb import make_lupton_rgb
from astropy.visualization import LogStretch
r1 = iband1[0].data
g1 = rband1[0].data
b1 = gband1[0].data
#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled1 = zscale(r1)
g_scaled1 = zscale(g1)
b_scaled1 = zscale(b1)
image = make_lupton_rgb(r_scaled1, g_scaled1, b_scaled1, Q = 1, stretch = .9, minimum = .3)
plt.imshow(image, origin = 'lower')
plt.show()
#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction
#based on whether this galaxy is spiral or elliptical.
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r1, (1050, 1050), (600, 600))
data0 = cutout.data
data = (data0-np.median(r1))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )
y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
p0 = [np.max(data), 10, 4, 300, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")
Fit: I_e = 1.13535, R_e = 17.60, n = 1.50, x0 = 266.66, y0 = 318.41
model_image = sersic_2d(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
iband1.close()
gband1.close()
rband1.close()
iband2 = fits.open("/Users/jackcolvin//m60_i-band_120.0s_bin1_250603_052218_jackcolvin_seo_0_FCAL.fits")
gband2 = fits.open("/Users/jackcolvin//m60_g-band_120.0s_bin1_250603_051942_jackcolvin_seo_0_FCAL.fits")
rband2 = fits.open("/Users/jackcolvin//m60_r-band_120.0s_bin1_250603_051649_jackcolvin_seo_0_FCAL.fits")
r2 = iband2[0].data
g2 = rband2[0].data
b2 = gband2[0].data
#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled2 = zscale(r2)
g_scaled2 = zscale(g2)
b_scaled2 = zscale(b2)
image = make_lupton_rgb(r_scaled2, .9*g_scaled2, b_scaled2, Q = 1, stretch = .9, minimum = .25)
plt.imshow(image, origin = 'lower')
plt.show()
#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction
#based on whether this galaxy is spiral or elliptical.
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r2, (1050, 1050), (600, 600))
data0 = cutout.data
data = (data0-np.median(r2))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )
y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
p0 = [np.max(data), 10, 4, 300, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")
Fit: I_e = 0.47759, R_e = 42.67, n = 2.00, x0 = 264.74, y0 = 325.45
model_image = sersic_2d(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
#again, it seems like the fit is not capturing the true steepness of the real curve
#in a very similar way to the first elliptical.
iband2.close()
gband2.close()
rband2.close()
iband3 = fits.open("/Users/jackcolvin//m49_i-band_120.0s_bin1_250603_051201_jackcolvin_seo_0_FCAL.fits")
gband3 = fits.open("/Users/jackcolvin//m49_g-band_120.0s_bin1_250603_050926_jackcolvin_seo_0_FCAL.fits")
rband3 = fits.open("/Users/jackcolvin//m49_r-band_120.0s_bin1_250603_050648_jackcolvin_seo_0_FCAL.fits")
r3 = iband3[0].data
g3 = rband3[0].data
b3 = gband3[0].data
#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled3 = zscale(r3)
g_scaled3 = zscale(g3)
b_scaled3 = zscale(b3)
image = make_lupton_rgb(r_scaled3, .93*g_scaled3, b_scaled3, Q = .5, stretch = .9, minimum = .4)
plt.imshow(image, origin = 'lower')
plt.show()
#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction
#based on whether this galaxy is spiral or elliptical.
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r3, (1050, 1050), (600, 600))
data0 = cutout.data
data = (data0-np.median(r3))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )
y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
y[i] += len(y[0]) - (i+1)
x[:,i] += i
p0 = [np.max(data), 10, 4, 300, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))
I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")
Fit: I_e = 0.34682, R_e = 54.06, n = 2.09, x0 = 255.71, y0 = 309.56
model_image = sersic_2d(x, y, *result.x)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)
plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()
#again, it seems like the fit is not capturing the true steepness of the real curve
#in a very similar way to the first elliptical.
iband3.close()
gband3.close()
rband3.close()
While I was unable to succesfully reach the first goal of this project – matching the expected sersic frofiles of n~1 for spiral galaxies and n~4 for ellipticals, I did still have some successes. After using a two component model for disk galaxies, I was able to measure a value of ~1 for 2/3 of my spiral galaxies, and n>2 for all three. Furthermore, I did measure elliptical galaxies to have a more dense light distribution than spirals, with all three eliptical galaxies having n~2.
I had to deal with significant sources of error throughput my analysis for both eliptical and spiral galaxies. In both cases, the tendency of SEO to spread light across a larger number of pixels due to its limited resolution (3-4 arcseconds) certainly had an impact on my fits. For spiral galaxies, the different light profiles of both the bulge and the spiral arms was a large source of errors that biased by values, and I was only able to partially mitigate those error source by a combination of masking, smoothing, and making a more complex, two-component model. For elliptical galaxies, The limited angular resolution of SEO likely had an even larger effect, as their light is much more similar to a point source, and thus requires high angular resolution to accurately measure. Furthermore, I observed a systematic underapproximation of sersic indeces in all three galaxies, evidenced by residual plots that sugested that a sharper model may actually be a better fit, at least for the central core of the galaxy.
As I continue to work on this project, I have a couple key considerations to keep in mind. First of all, I want to continue tweaking my two-component spiral galaxy model to try to retrieve a model that can simulataneously give me accurate values for both my bulge and my disk. Furthermore, I wasnt to work more on smoothing to lessen the effect that the spiral arms may be having on my sersic fits. For the eliptical galaxies, I want to try to contrain my other parameters apart from sersic index more tightly to see if i can retrieve a model that qualitatively reduced my residuals more effectively. Furthermore, for both cases, I want to both take more SEO eposures, and maybe also try some higher reolution images from the Legacy or SDSS surveys, so try to isolate whether my errors are a result of SEO's poor resolution, or some aspect of my specific analytical techniques.